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(Abstract) In this paper we extend the cubic B-splines collocation method to enable it to solve one-dimensional hyperbolic 
(eigenvalue and wave propagation) problems under arbitrary boundary conditions. This is achieved by analogy with the finite 
element method, introducing a 'collocation mass matrix' that cooperates with the previously known system matrix, now called 
'collocation stiffness matrix'. In agreement with earlier findings on elliptic problems, we found that in time-dependent 
problems it is again sufficient to use double internal knots (C'-continuity) in conjunction with two collocation points between 
successive breakpoints. In this way, the number of unknowns becomes equal to the number of equations, which is twice the 
number of breakpoints. We paid particular attention to the handling of Neumann-type boundary conditions, where we found it 
necessary to properly eliminate a column in both mass and stiffness matrices. For the first time, we found that the cubic B-splines 
collocation procedure (with C 1 -continuity) leads to identical results with those obtained using piecewise Hermite collocation. The 
numerical examples show an excellent quality of the numerical solution, which is far superior to that of the conventional finite 
element method, for the same number of nodal points. 
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1. INTRODUCTION 

In the framework of global interpolation methods based on 
Lagrange polynomials, it has been recently shown that the 
computation effort may be substantially reduced when 
applying the collocation instead of the Galerkin-Ritz 
formulation [1]. However, using spline curves, or piecewise 
polynomials, is more effective in representing the solution to 
the differential equation than pure polynomials [2]. The 
volume of Ascher et al. [3] provides a treatise on spline 
bases, collocation theory, and spline collocation for 
application to the numerical solution of 
boundary-value-problem (BVP) for ordinary differential 
equations (ODE). Fairweather and Meade [4] give an 
extensive review (273 papers covering the period 1934-1989) 
of collocation methods and various implementations. They 
describe the most common forms of collocation, including 
nodal, orthogonal, and collocation/Galerkin. An early work 
towards the solution of eigenvalue problems, however based 
on Schoenberg's formulation [5], is [6]. 

In spite of the abovementioned work, so far most part of the 
relevant research focuses on mathematical topics and reduces 
mainly to elliptic problems. For example, the FORTRAN 
codes of [2] have been implemented also in the MATLAB 
software (formerly under the name splines tool) [7], but they 
operate 'as is' for the solution of linear and nonlinear elliptic 
problems only. This happens because (i) the eigenvalue and 
transient analysis require the construction of mass and 
stiffness collocation matrices, and (ii) although the 



implementation of Dirichlet type boundary conditions is a 
trivial task, the same does not hold for Neumann-type ones 
that require a special treatment. As for the state-of-the-art, the 
implementation of the function 'spcol' in time-dependent 
analysis has been made in very few biometrics [8] and 
chemical engineering applications [9]. 

In this context, the primary aim of this paper is to 
investigate the applicability and performance of deBoor's 
methodology [2] in the numerical solution of hyperbolic 
problems (eigenvalue and transient). The main novel feature 
of this work is the development of mass and stiffness 
collocation matrices analogous to the finite element method. 
Consequently, standard eigenvalue analysis based on the QR 
algorithm, and standard time-integration techniques such as 
the central difference method will be applied. 

In addition to the abovementioned global cubic B-splines 
collocation, piecewise-Heimite polynomials (without upwind 
features) that act between adjacent breakpoints [10] will be 
compared for the first time. The finding of coincidence is the 
secondary novel feature of this work. 

The theory is sustained by four one-dimensional numerical 
examples from the field of applied mechanics. 

2. FORMULATION 

Below is the formulation of typical static (thermal) and 
dynamic problems and then follows the global B-spline and 
piecewise Hermite interpolation. 

2.1 Thermal Analysis 
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The steady-state thermal behavior of structures follows 
Laplace equation. In the particular case of a radially 
symmetric structure the PDE degenerates to a simple ODE 
in terms of the radius r as follows: 



d u 1 du 
— + — = 

dr r dr 



(1) 



2.2 Elastodynamics - Wave equation 

The vibration of an elastic rod is described through the usual 
hyperbolic differential equation: 



d_ 
dx 



du 
EA — 

V dx 



d 2 u 
dS 



(2) 



where E and p denote the material properties, i.e. the 
Young's modulus and the mass density, respectively, and A 
denotes the cross sectional area of the rod. Finally, / 
denotes the distributed load along the longitudinal x-axis, 
while t denotes the time. In dynamic analysis, the axial 
displacement, u, is a function of both space and time (i.e. 
u — u(x,t\), while in static analysis it is only a function of 
space (i.e. w=w(x)). 

Introducing the velocity of the longitudinal elastic wave: 



Je/p> 



(3) 



in the absence of body forces ( / = ), Eq.2 receives the 
well known form of an acoustic wave: 



(l/c 2 )-d 2 u/dt 2 -d 2 u/dx 2 = 0. 



(4) 



It is well known that, besides elastic waves, Eq.4 can also 
describe wave propagation within acoustic pipes of constant 
cross-sectional area; in such a case the variable u corresponds 
to the acoustic pressure. 

2.3 B -splines Interpolation 

B-splines interpolation can be found in many textbooks 
[2,11], in which efficient procedures to determine the 
involved basis functions, N. \x) , can be found. In general, 
the interpolation within an interval < X < L is given by: 



(5) 



Concerning Eq.5, it should become clear that: 
The starting point for the construction of the functional 
set, N j , is the selection of the breakpoints and the 
piecewise polynomial degree, p. Let us assume m 
breakpoints, which obviously define (m-1) segments. In the 
sequence, let the polynomial degree be p = 3 (cubic 
B-splines). 

(i) The variables a. involved in Eq.5 are not directly 
associated to nodal values (those at the breakpoints), 



except of those two corresponding to the end points. 

(ii) The value of 'n' depends on the multiplicity chosen for 
the internal nodes. In this paper, we require C l -continuity, 
thus the multiplicity must be equal to two; therefore, 'n' is 
twice the number of breakpoints (n = 2m). In general, the 
spline curve is C p ~ -continuous everywhere except at the 
location of the repeated knots where it is 
C F ~ -continuous. 

(iii) Therefore, if we choose two collocation points within 
any of the aforementioned (m- 1 ) segments defined by the 
breakpoints, we can obtain so many equations as the 
number of the unknown coefficients. In case of elliptic 
problems, the best choice was found to be two Gauss 
points at the well-known relative positions ( g = ±1/ \3 ), 
as also was found by deBoor and Swartz [12] (in elliptic 
problems). 

Typical plots of basis-functions are shown in Figure 1, in 
which double knots have been considered. It is noted that both 
the first and the last basis functions correspond to the left-end 
and the right-end axial displacement, respectively; these are 
also a mirror of each other with respect to the center of the 
interval [0,L], and equal to unity at the two ends. 




□ 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.S 0.3 1 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 1. Basis functions in the domain [0,1]. In the top and the 
bottom there are four subintervals (m = 5) and eight subintervals (m = 
9), respectively; in both cases the internal breakpoints are considered 
with multiplicity of internal knots equal to two. 

2.4 Piecewise-Hermite Interpolation 

The domain is again uniformly divided into (m-1) segments 
using m nodal points (the ends are included). Starting from left 
to the right, the nodes are numbered by ascending order, 
whereas the 2m degrees of freedom are as follows: 
(u lt q \(u ,q ~),...,(u ,q ") with g = (Sw/dx) . In this 
way, having considered the flux as an independent DOF, 
C'-continuity is ensured. 

Between any of the (m-1) segments, in local numbering the 
variable u is approximated by: 
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u(x,t) = AT (x)u i (t) + N 2 (x)q l (t) 
+ N 3 (x)u 2 (t) + N i (x)q 2 (t) 



(6) 



where the shape functions are the well-known Hermite 
polynomials (Figure 2), which are expressed in terms of the 
normalized coordinate w = (x/L) '■ 



AT ( w ) = l-3w 2 +2w 3 =1-3 
jV 2 (w) = w — 2w 2 + w 3 = 
A^ 3 (w) = 3vv 2 - 2w' = 3 
N 4 (w) = -w" + w 3 = 



f- 



f- 



- -2 - I - 



A" 



u 



f- 



(7) 













kl) _ 



It is noted that the above notation presumes that the fluxes q 
and are both positive when directed to the positive x-axis. 
Then q causes 'compression' while q 2 causes 'tension'; 
also the corresponding shape functions in Figure 2 become 
anti-symmetric with respect to the middle of the interval [0,L] 
and obtain positive and negative values, respectively. If the 
convention for q 2 changes, both shape functions may 
become positive and symmetric with respect to the middle of 
the domain. Then, the 'mirror' property is obtained for both 
translational ( u x ,u ) and rotational ( q { ,q 2 ) degrees of 
freedom. 




D D.I 0.2 D.3 4 0.5 06 7 08 09 1 



Figure 2. Hermite polynomials defined in [0,1] (blue: translational 
DOF (N 1 jy 3 \ green: rotational DOF (N 2 ,N 4 )). 

3. COMPUTATIONAL PROCEDURE 

In all cases below, the ODE is collocated at the n = 2m Gauss 
points in the interior of the subintervals, i.e. at the locations: 

(x.+x.) 1 (x, -x. . ) , \ 

2 V3 2 

Below, the computational procedure for dealing with 



eigenvalue and wave propagation problems is presented. 
3.1 Eigenvalue Problem 

Collocating Eq.4 at the abovementioned n = 2m Gauss points 
one obtains the well known matrix formulation [1]: 



[M]{a(0} + [K]{a(0} = {f(0} 



(9) 



where [M] and [K] , both of dimensions 
(«-2)xhe 2[m — l) x 2m , are the nonsymmetric mass and 
stiffness collocation matrices, respectively, which are given 
by: 



r l ^ 



m — 

•j 



2 



(10) 



in 



fc..=-Ar;(x), i = l,...,2(m-l), j = l,...,2 

with (x) denoting the shape functions (in global 
numbering) at the collocation point x. , in ascending order 
from left to the right. 

It is remarkable that Eq.9 covers both cases, i.e. the global 
B-splines interpolation (cf. Eq.5) as well as the piecewise 
Hermite interpolation (cf. Eq.6 and Eq.7). 

Concerning the eigenvalue (free vibration) problem, the 
force vector |f (f )| is taken equal to zero and, therefore, the 
eigenvalues are extracted by requiring: 



det||[K]-« 2 [M]| = 



(11) 



In this study we used the standard MATLAB function 'eig'. 
3.2 Wave Propagation Problem 

The wave propagation problem can be generally solved using 
either modal analysis or one of the well known time-integration 
methods (explicit, implicit). In this paper we choose the central 
difference method [13]: 

1 , , 1 , , 
u = lu — 2u +u j ) and u = \u — u ). (12) 

At 2 " " 2At 

Substituting Eq.12 into Eq.9, the last written for 
t , — t + At , one obtains: 



f 1 > 

M 

vAf 2 j 



K 



M 

At 2 j 



f 1 > 
M 

\At 2 j 



(13) 



from which we can solve for u . . 

n + l 

3.3 Boundary Conditions 

Despite the apparent simplicity of the procedure, some 
technical difficulties appear and require special treatment and 
explanation. First of all, it is reminded that two DOF exist at 
every nodal point (or breakpoint) for the piecewise cubic 
Hermite (or B-splines) collocation formulation. Concerning 
the left and right ends, the DOF correspond to the quantities u 
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and q = ( dujdx ) m . In every formulation, the vector of the 
coefficients is as follows (it is reminded that the domain is 
uniformly divided into (m-1) segments using m nodal points, 
where the ends are included): 
Piecewise cubic Hermite: 

{u l ,q r u 2 ,q 2 ,...,u m _ l ,q m _ l ,u m ,qJ T 
Cubic B-splines: 

= Mj, Cl 2 , a 3 , Cl 4 , . . . , <Z 2m _ 3 , 2m -2 > fl 2m-l > fl 2m = M m } 

Therefore, for a bar fixed at both ends (or an acoustic pipe 
open at its both ends), in the piecewise Hermite formulation it 
is sufficient to eliminate the first ( U x = ) and the last minus 
one (u = 0) columns of the matrices [M] and [K], while in 
B-splines formulation we must eliminate the first and the last 
column (it is reminded that B-splines have the "mirror" 
property). 

The situation becomes more complex when dealing with a 
free end (Neumann-type), for example at x=L. In this case the 
strain (or flux, or velocity) is zero or has a given value, but it is 
not always allowed to eliminate the corresponding column. In 
more details: 

(i) In piecewise Hermite formulation the DOF that 
corresponds to the vanishing flux at the end is quite 
independent of the rest DOFs, and therefore it can be 
directly eliminated. In contrast to the FEM analysis, none 
row is eliminated but only the corresponding columns. 
Therefore, starting from a matrix system of 
(n - 2) x n = l{m - 1) x 2m dimensions, the elimination 
of two DOF leads to a system of 
(n - 2)x (n - 2) = l{m -l) x 2(m-l) dimensions. 

(ii) In B-splines formulation the elastic strain (or the acoustic 
flux/velocity) at the right free end (x=L) is not an 
independent variable but it can be derived by taking the 
first derivative of Eq.5: 

4(L) = M '( L ) = Z< P ( L )- fl , (14a) 

1=1 

Due to the compact support of the basis functions, for p=3 
only the last two out of the n=2m basis functions have a 
non-vanishing derivative at x=L. Therefore, it holds: 

q(L) = N' n _ l {L)-a n _ l +N' n {L)-a n 

[q(L)-N'(L)-a] (14b) 
=> a = 

KM 

For the solution of the eigenvalue problem we can assume 
that g(L) = , and therefore Eq.l4b is somehow 
simplified. Concerning the matrix term 
([k] - a> [m]) • {a} , every row can be transformed from 
its initial form: 
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\k a + ... + k. a +k. a +k. a) 

... . . ,.. . (15) 

-co' (m n a l + ... + m. n _ 2 a u 2 + m. u _fl n l + m. a n ) 
to the final expression: 

[ka + ... + k,_a n _ 2 + (* „ - AT (L)/A^ (L) • k ) « ] 
-co 2 [m ti a +... + m. n t a^ + (m. n -N' n (L)j N' n { (L) • m. n _, )aj 

(16) 

Therefore, concerning again B-splines collocation, in case of 
one Dirichlet boundary condition at x=0 and one Neumann 
condition atx=L, we proceed as follows. First we eliminate the 
first column ( a y —u l — ). Then we condense the last two 
columns in one by subtracting the (n-l)-th column (multiplied 
by the coefficient N' n (L) / N'^L)) from the last one. In 
this way we finally derive a system of 
(ti-2)x(ti-2) = 2(m-l)x2(m-l) dimensions. 

Note that the abovementioned condensation is not 
applicable to the piecewise Hermite formulation, obviously 
because N' 3 (l) = 0, N[(w) = l. 

4. NUMERICAL RESULTS 

The performance of the proposed global collocation method 
will be evaluated in four characteristic test cases: one static 
(to show the convergence quality) and three dynamic ones. As 
previously mentioned, from an engineering point of view the 
dynamic problems may concern either an elastic rod of 
constant cross-sectional area or, equivalently, an acoustic 
straight pipe. Since previous studies [12,14] concerning the 
elliptic problem have shown that the best choice is to use 
cubic B-splines in conjunction with collocation points taken at 
the two Gauss points between the uniformly distributed 
breakpoints, results will be presented for these conditions 
only. For comparison purposes, the collocation method is 
compared with the conventional finite element (FEM) 
solution, using the same nodes as the breakpoints. 

4.1. Problem 1: Steady-state Conduction in a 
Cylindrical Wall 

This example was chosen because it possesses a 
non-polynomial (logarithmic) solution. It refers to the 
potential problem (Laplace equation) with a "steep" 
temperature distribution, occurring along the radius of a thick 
long hollow cylinder. The latter is subject to given uniform 
inner surface temperature (u { = 1000 °C) and a given uniform 
outer surface temperature (u 2 = °C) while the inner and outer 
radii are Ri = lm and R 2 = 32m, respectively. The analytical 
temperature distribution depends on the radial direction only: 

u (r) = «, + (u 2 - tij )/ln (tf 2 /*, ) • In ) ( 17 > 

The radius is divided into a variable number of n = 4, 8, 16, 
32 and 64 uniform segments. Based on the normalized L 2 
error norm ( L u = L, x 100 in %): 
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The excellent quality of/the convergence is shown in Table 1. 
Table 1. Convergence quality of the error norm. 



Number of subdivisions 

(») 


Energy norm 

L u (in %) 


4 


1.19E00 


8 


0.73E-01 


16 


0.20E-02 


32 


2.54E-05 


64 


1.74E-07 



4.2 Problem 2: Elastic Rod Fixed at Both Ends 

A rod of length L is fixed at x = and x = L and is subject to 
axial free vibration (equivalently, we could refer to an acoustic 
pipe with both ends open). The exact eigenvalues are given by: 



co] =(ix/L) 2 -(E/p), i = l,2,. 



(19) 



The rod is uniformly divided into 1, 2, 4, 8, 16, 32 and 64 
segments. The numerical results are shown in Figure 3 where 
one can notice that the convergence is monotonic, of high 
quality and far superior to the conventional FEM solution. 
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Figure 3. Errors (in %) of the calculated eigenvalues for an elastic rod 
with fixed ends (the horizontal axis refers to the number n of elements 
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Figure 4. Errors (in %) of the calculated eigenvalues for an elastic rod 
with fixed left and free right end (the horizontal axis refers to the 
number n of elements or subintervals of breakpoints). 

4.3 Problem 3: Elastic Rod Fixed at Left and Free 
at Right End 

The same rod is now analyzed for different boundary 
conditions (the equivalent acoustic pipe is open at left and 
close at right end). The quality of the calculated eigenvalues 
and the superiority to the FEM is shown in Figure 4, where it is 
expressed as an error (per cent) with respect to the exact 
solution: 



co*=[(2i-l)7r/(2L)f-E/p, 1 = 1,2,. 



(20) 



4.4 Problem 4: Elastic Rod Subjected to a 
Heaviside Load 

An elastic rod of length L is fixed at one of its extremities 
(x=0) and it is subjected to an axial Heaviside type loading 
a L = Eq Q H (t-0) [N/m 2 ] at the other one (x=L). The 
analytical solution can be found in [15,16]. For simplicity, all 
geometric and material data were assigned the unitary value. 
The (explicit) central difference scheme was appl ied. 

For the case of n=8 uniform intervals ( c = ^Je/ p =lm/s), 
where the distance between two successive nodes equals to 
Ax= 0.125m, a very upper limit for the time step could be 
At a =0.125 (CFL-criterion). However, for the sake of 
conservatism, the time step has been conservatively chosen 
equal to At = 0.02 Ax/c =0.0025 s. The time history (6000 
steps) of the normalized axial displacement at the free end 
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(x=L) and the middle point (x=L/2) is shown in Figure 5 and 
is of excellent quality. 



2 



1.5 







-05 




ct/L 



Figure 5. Axial displacement at the free end (x=L) [blue line] and the 
middle [red line], using the proposed collocation method n = 8 
uniform intervals and time step At = 0.2 Ax/ c = 0.0025 s 

5. DISCUSSION 

Similar to previous findings that were obtained for 
polynomials [1], the superiority of the proposed collocation 
method to that of the finite element method is due to its global 
character and the relatively high degree of the piecewise 
polynomials (p = 3). Clearly, the cubic B-splines interpolation 
is very smooth and can accurately approximate the involved 
eigenvectors, while the conventional finite element method 
uses piecewise linear (hat) shape functions that appear sharp 
edges at their interfaces (they posses only C°-continuity). 

In all elliptic and hyperbolic numerical examples of this 
study we found that the cubic B-splines collocation method 
leads to identical results with those obtained using the 
piecewise cubic Hermite collocation method. This numerical 
finding can be theoretically explained as follows. Since two 
knots are associated to each breakpoint (multiplicity equal to 
two) it implies that the splines curve is generally 
C p 2 -continuous, where p is the degree of the piecewise 
polynomial. But as in this study the polynomial degree was 
chosen to be cubic ( p — 3 ), the splines curve is 
C 1 -continuous, a fact that is consistent with the 
C 1 -continuity involved in piecewise Hermite approximation 
(also of third degree). 

Although in this type of problems the proposed B-splines 
interpolation is identical (as for the quality of the obtained 
results) with the piecewise cubic Hermite collocation method, 
the first is the vehicle to choose any multiplicity for the 
internal knots. For example, although no relevant results were 
reported in this study, we could alternatively choose the 
multiplicity be equal to one, instead of two. In such a case the 
numerical solution would be C 2 -continuous, as it would 
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include a full cubic polynomial plus some Schoenberg's 
truncated cubic monomials [5]. For example, if we divide the 
domain [0,L] into four segments, cubic B-splines 
approximation would lead to seven control points and seven 
associated coefficients, a., i = 1, . . . , 7 . However, since we 
deal with a two-point BVP, we have to determine only five 
out of the seven aforementioned coefficients. To this purpose 
one must choose at least five collocation points within the 
domain of length L. Obviously, one could use even more than 
five collocation points but then, since the number of equations 
would be larger than the number of variables, a least-squares 
scheme had to be applied. In the same example, although a 
choice of five (or more, in combination with the least-squares 
method) Gauss points or Chebyshev roots generally works, 
the quality of results becomes questionable while the 
DeBoor's approach (double knots, and two collocation points 
between successive breakpoints) is always mathematically 
robust [12]. 

6. CONCLUSIONS 

In this study we achieved to extend the well-known cubic 
B-splines collocation method, previously used in the solution 
of elliptic problems only, to one-dimensional hyperbolic ones 
under arbitrary (Dirichlet of Neumann) boundary conditions. 
The key point was that, in addition to the system matrix, we 
achieved to create a mass collocation matrix; thus the 
computational problem was formulated very similarly to the 
well known finite element method. Consequently, the 
eigenvalues were calculated in an algebraic way using standard 
QR algorithms, while wave propagation (time response) was 
predicted using a simple central-difference scheme for 
time-integration. The findings of this study show that the 
proposed technique has an excellent convergence and its 
overall performance is far superior to the conventional finite 
element method for the same number of nodal points. Another 
interesting finding is that, B-splines collocation based on two 
knots per breakpoint leads to identical eigenvalues and time 
response with those obtained using again the collocation 
method but in conjunction with piecewise cubic Hermitian 
polynomials. Although the results of this study reduced to ID 
problems only, on-going research has revealed that extension 
to 2D (quadrilaterals) and 3D (boxlike) time-dependent 
problems is a straightforward procedure on the basis of tensor 
products of the involved control points per direction. 
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